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ABSTRACT 

Results are presented from a time-dependent, numerical investigation of spherical 
accretion onto black holes, within the framework of relativistic radiation hydrody- 
namics. We have studied the stability of self-consistent, stationary solutions of black 
hole accretion with respect to thermal and radiative perturbations and also the non- 
linear evolution of unstable, high temperature models, heated by the hard radiation 
produced by the accretion flow itself in the inner region near to the horizon. In some 
cases, a hydrodynamic shock forms at around 10 3 -10 4 Schwarzschild radii, where 
Compton heating exceeds radiative cooling. The calculations were made using a suit- 
ably designed radiation hydrodynamics code, in which radiative transfer is handled 
by means of the PSTF moment formalism and which contains an original treatment 
of the radiation temperature equation. 
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1 INTRODUCTION 

Stationary, spherical accretion onto black holes is a well- 
known and extensively studied topic. Starting from the sem- 
inal paper by Bondi (1952), many papers have been devoted 
to the analysis of spherical accretion under a variety of con- 
ditions, mainly in order to obtain a definite estimate of the 
efficiency of the process. In contrast with accretion onto neu- 
tron stars, the efficiency is not fixed by the requirement that 
all the kinetic energy of the accreting gas must be converted 
into radiation, since no rigid surface exists which can stop 
the flow. Matter can cross the horizon carrying a substan- 
tial fraction of the gravitational potential energy liberated 
and the efficiency of the process is determined solely by the 
effectiveness of the radiative processes in converting the in- 
ternal energy of the accreting gas into radiation, as first 
noted by Shvartsman (1971). As shown by many authors 
(see e.g. Michel 1972; Novikov & Thorne 1973; Blumenthal 
& Mathews 1976; Begelman 1978; Brinkmann 1980), the 
flow properties are fixed once the accretion rate is specified, 
so that stationary solutions can be completely characterized 
by their position in the (rh, I) plane, where rh and I are, 
respectively, the accretion rate and luminosity in Edding- 
ton units. Stationary spherical accretion onto black holes 
was investigated in detail by Nobili, Turolla & Zampieri 
(1991, hereafter NTZ); Fig. 1 shows the (rh, I) diagram for 
the complete set of their solutions. At low accretion rates 
(rh < 1), spherical accretion is very inefficient in convert- 
ing gravitational energy into radiation since the density is 



too low for cooling processes to be effective; the emitted lu- 
minosity is also very low (lower branch in Fig. 1, hereafter 
the LL branch). These models are essentially adiabatic and 
have very high temperatures (see also Shapiro 1973a). In 
this regime, the only possibility for increasing the efficiency 
of the accretion process is related to the presence of mag- 
netic fields, which can cause strong dissipation (e.g. through 
reconnection of field lines) and induce strong emission of 
synchrotron radiation (Shapiro 1973b; Meszaros 1975). Sof- 
fel (1982) studied in some detail the transition from the 
optically thin regime to the optically thick one: as rh in- 
creases, the cooling processes become more effective and the 
gas temperature decreases, causing in turn a decrease in the 
total emitted luminosity (with a local minimum at around 
rh ~ 0.1). For higher values of the accretion rate, free-free 
absorption is no longer negligible and the gas becomes op- 
tically thick in the inner region near to the horizon of the 
black hole. The temperature increases because heating ex- 
ceeds cooling and also the luminosity rises since radiation 
is in LTE with the gas in the inner core. Preliminary inves- 
tigations of spherical accretion in the diffusion regime were 
made by Tamazawa et al. (1974), Maraschi, Reina & Treves 
(1974), Kafka & Meszaros (1976), Vitello (1978), Begelman 
(1979), Gillman & Stellingwerf (1980) and Freihoffer (1981), 
while a complete treatment was finally given by Flammang 
(1982, 1984), who showed the existence of a subcritical point 
related to the equation for the radiative luminosity. When 
rh > 1, the inner core starts to be optically thick to electron 
scattering as well. In these conditions, a trapping radius ap- 
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pears (Rees 1978; Begelman 1978), below which photons are 
advected into the black hole, since their outward diffusion 
velocity is smaller than the inward velocity of the accretion 
flow. This makes the process less efficient and the rate of 
increase of luminosity with rn becomes slower. The station- 
ary, hypercritical regime has been investigated by Blondin 
(1986). 

For 3 Ss m Ss 100, there is also another class of solu- 
tions, characterized by having high temperatures and lumi- 
nosities (upper branch in Fig. 1, hereafter the HL branch). 
These are dominated by the effects of comptonization which 
keeps the gas and radiation temperatures almost equal in 
the inner part of the flow where the density is sufficiently 
high to make the Compton parameter very large. In the in- 
termediate region between 10 2 and 10 5 r g (where r g is the 
Schwarzschild radius of the black hole), Compton heating 
dominates and the only competitive cooling mechanism is 
free-free emission. The first authors to investigate the pos- 
sible existence of high luminosity solutions were Wandel, 
Yahil & Milgrom (1984) and Park (1990a,b), who performed 
a detailed study of spherical accretion for a large range of ac- 
cretion rates and considered also a two-temperature model 
with pair production. High luminosity stationary solutions 
have relatively high efficiency and appear to exist only for a 
very definite range of accretion rates. Already in 1976, Os- 
triker and collaborators (Ostriker et al. 1976) pointed out 
that, because of the non-local nature of comptonization, the 
heating produced in the flow by the hard radiation coming 
from the inner region can increase the gas temperature in 
such a way that the internal energy density becomes larger 
than the gravitational energy density and the accretion pro- 
cess is then stopped. This effect is called preheating. Later 
on, Cowie et al. (1978), Shull (1979), Stellingwerf & Buff 
(1982) and Krolik & London (1983) showed that preheating 
is very important in placing limits on the region of parame- 
ter space within which HL solutions for black hole accretion 
can exist, although the strength of preheating is reduced if 
Compton cooling is taken into account (see e.g. Bisnovatyi- 
Kogan & Blinnikov 1980). As shown by NTZ, preheating at 
the sonic point for the matter prevents the existence of high 
luminosity solutions with m <, 3, while preheating within 
the sonic radius prevents the existence of stationary solu- 
tions for m J> 100. 

The stability of these solutions remains a completely 
open question and the main goal of the present paper is 
to study this. The first attempt to investigate the stabil- 
ity of isothermal accretion was made by Stellingwerf & Buff 
(1978) using an Eulerian scheme, based on an extension of 
the Henyey relaxation method. They found that transonic 
accretion is quite stable. By means of a general relativis- 
tic analytical calculation, Moncrief (1980) showed that, for 
isentropic flows, no unstable normal modes exist which ex- 
tend outside the sound horizon. The first studies includ- 
ing the heating and cooling terms due to the presence of 
the radiation field (Cowie et al. 1978; Stellingwerf & Buff 
1982; Stellingwerf 1982) were devoted to analysing the ef- 
fect of preheating on the stability of the accretion flow and 
to defining the region of the (rh, I) plane where the exis- 
tence of stationary solutions is not allowed. In particular, 
Stellingwerf (1982) presented a local stability analysis of op- 
tically thin, X-ray heated accretion flows and showed that, 
for sufficiently high luminosities, a finite amplitude drift in- 



stability can develop, due to the form of the free-free cooling 
function, causing a time-dependent behaviour of the solu- 
tion on a time-scale ranging between a day and a few tens 
of days. Krolik & London (1983) used the WKB method 
to derive the dispersion relation for modes coupling den- 
sity, temperature and velocity perturbations in an optically 
thin, accreting gas and found that, although stationary solu- 
tions with high temperature and luminosity can exist, heat- 
ing of the gas inside the sonic radius leads to the onset of 
a thermal instability in a large region of the (m, /) plane. 
Gilden & Wheeler (1980) and Vitello (1984) investigated 
time-dependent, optically thick accretion within the frame- 
work of General Relativity, treating the radiation field in 
the diffusion approximation and using two different numer- 
ical schemes: a Lagrangian hydrodynamic code in the first 
case and a Linearized Block Implicit Algorithm in the second 
one. They found that, within this approximation, no matter 
which initial conditions the code started from, convergence 
toward stationary LL solutions was rapidly achieved, show- 
ing that they are intrinsically stable. Finally, NTZ, using an 
argument originally suggested by Nobili, Calvani & Turolla 
(1985) and based on Prigogine's criterion, argued that HL 
solutions might be unstable because of the large value of the 
entropy production rate. 

Despite the fact that the stationary problem has been 
extensively investigated, mainly for shedding light on the 
efficiency of the radiation generation, we think that several 
aspects require further consideration such as, for instance, 
investigating the stability properties of the high luminosity 
solutions and searching for the existence and non-linear evo- 
lution of possible heated or shocked models. In the present 
paper we present an analysis of the stability and time- 
dependent behaviour of the solutions for spherical accretion 
onto black holes within the framework of general relativis- 
ts radiation hydrodynamics. We adopt the conventions that 
Greek indices run from to 3 and covariant derivatives are 
denoted with a semi-colon; a spacelike signature ( — ,+,+,+) 
is used. 



2 THEORY 

In this section we review the derivation of the equations 
of relativistic radiation hydrodynamics in spherical symme- 
try for a self-gravitating matter fluid which is interacting 
with radiation. A complete treatment of this subject was 
recently presented by Rezzolla & Miller (1994) in connec- 
tion with the cosmological quark-hadron transition (see also 
Schmid-Burgk 1978, Mihalas & Mihalas 1984 and Park 1993 
for a discussion of the equations of radiation hydrodynam- 
ics in spherical flows) . Since the source terms which account 
for the emission and absorption of photons are more easily 
written in the reference frame comoving with the gas, this 
will be our fiducial frame and u a will denote its 4-velocity. 
In this frame the energy and momentum conservation equa- 
tions can be obtained by projecting the 4-divergence of the 
stress-energy tensor of the matter plus radiation fluid along 
u a and orthogonal to it, giving 

«a(T^+T^) ;/3 =0, (1) 

M T m +r;f) ;/3 = o, (2) 
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Figure 1 
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Figure 1. The (m, I) diagram for the complete set of stationary solutions found by NTZ (circles). Also shown are the 6 initial models 
whose relevant parameters are listed in table 1. Filled circles mark the stable LL stationary solutions, while open circles denote the 
unstable HL models. Asterisks indicate the low fa HL solutions, which might still be unstable, but on much longer time-scales. 



where h la — g la +UjU a is the projection tensor orthogonal 
to the 4-velocity and 

T^^(e+p)u a ^+pg^, 



1 



(3) 



(4) 



are the matter and radiation stress-energy tensors, respec- 
tively; e and p are the energy density and pressure of the 
gas. This system of equations must be supplemented with 
the rest-mass conservation equation 



(5) 



where p is the rest-mass density of the matter measured in 
the comoving frame. In equation (4) the stress-energy tensor 
of the radiation field is expressed in terms of the Projected 
Symmetric Trace-Free (PSTF) moments (Thorne 1981) 



In ai ...n ak dQ 



Trace — Free 



(6) 



where / = I(x a ,p a ) is the specific intensity of the radia- 
tion field and n a is the unit vector which gives the direc- 



tion of propagation of a photon as seen in the rest frame 
of the fiducial observer 11°. Integration is over solid angle 
in the projected space and "Trace-Free" denotes the con- 
sequence of the usual tensor operation. By definition, the 
PSTF moments are symmetric tensors which lie entirely in 
the projected space and represent the relativistic analogue 
of the classical moments of the specific intensity (see e.g. 
Chandrasekhar 1960). In terms of PSTF moments, / can be 
written as 



fc=0 



(2fc + l)!! 
47rfc! 



M c 



(7) 



The specific intensity obeys the general relativistic equation 
of radiative transfer 

^ = s - w 

where iV = (c 2 /2h)(I /v^) is the photon occupation num- 
ber (y = cu a p a /h is the photon frequency measured in the 
comoving frame) , I is a non-affine parameter along the pho- 
ton trajectory in phase-space and S is a source function 
which describes the effects of the interaction between mat- 
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ter and radiation, its actual form depending on the radiative 
processes which are considered. The moments of the source 
function, S ai "' ak , can be defined in analogy with equation 
(6). If there is no interaction, S — and N is conserved along 
each photon trajectory. Moment equations can be obtained 
from equation (8) by inserting the expansion (7) for I (and 
the equivalent one for S) and taking the PSTF part (i.e. pro- 
jecting orthogonal to u a , removing the trace and performing 
the symmctrization) ; this gives rise to an infinite hierarchy 
of differential equations. 

Finally, in order to calculate self-consistently the metric 
tensor g al3 which describes the geometry of the space-time, 
we need to solve the Einstein Field Equations for the system 



R° 



-~2 9 



a/3 



R - 



BttG 



\ M + 1 R 



(9) 



In spherical symmetry we can define a local orthonor- 
mal frame comoving with the flow as {eg, , e^, e^}, with 
e ( -, = u, ef- being in the radial direction and eg, being 
orthogonal to each other and to . Since in this case / and 
S are invariant under rotations of the photon direction n a 
about ef, it is possible to show that all of the components 
of each PSTF moment of rank k can be evaluated in the 
comoving frame as functions of the radial one 



w k = M r 



Sk 



2tt 



k\(2 k 
~(2k 



2tt 



k\(2 k 
(2k 



fc + l)l frpk 

+ 1)!! cj 

k + l) h f , 
+ 1)!! c 3 J 



(n)dp, 



SP k {p)dp, 



(10a) 



(106) 



where P k {p) is the Legendre polynomial of order k, p = n r 
and a denotes the e a component. In particular 



M = Wo .. 



M a = wie?, 



M 



a/3 ( 
= W2 I 



a B 



a 13 
-C; e- 
2 e 



1 a 3 

2 <p V 







(11) 

(12) 



(13) 



In the following, we will be interested only in frequency- 
integrated moments and we will use wu to denote the inte- 
gral over v of w k (v). It is easy to see that wo is the radiation 
energy density and cwi is the radial component of the ra- 
diative flux. 

We next introduce the spherically symmetric, comoving- 
frame line element 



ds 2 = —a 2 c 2 dt 2 + b 2 d[i 2 + r 2 (d6 2 + sin 2 6dip 2 ^ , 



(14) 



where t and n are the Lagrangian time and the comoving ra- 
dial coordinate (taken to be the rest mass contained within 
a comoving spherical shell), r is the Eulerian radial coordi- 
nate and a and b are two functions of t and /i which need 
to be computed. The complete system of radiation hydro- 
dynamic equations (1), (2) and (5) along with the Einstein 
Field Equations (9) can be cast in the form (Rezzolla & 
Miller 1994) 



et — hpt + acso = Energy equation , 



u t + ac 
GM 

' 9 9" 



£ / p M + bsi \ 4nGr 
b I ph ) c 4 



(p + + 11)2^ 



(15) 



(16) 



= Euler equation , 



pr 2 



b = 



1 Un — 4irGbrwi lc 
+ ac' 



= Continuity equation , 



1 



Airr 2 p ' 



(a/i) M hpp - e M + bsi 



ah 



hp 



= 0, 



where 

n 

u = — 

ac 



(17) 



(18) 



(19) 



(20) 



(21) 



is the radial component of the fluid 4-velocity measured in 
the fixed Eulerian frame, V = (1 + u 2 -2GM/c 2 r) 1/2 = r M /6 
is the general relativistic analogue of the Lorentz factor, M 
represents the effective gravitational mass (for black hole + 
gas + radiation) contained within radius r and h = (e + 
V)l P is the specific enthalpy. Here, the subscripts t and p 
denote partial derivatives with respect to the corresponding 
variables and so and si are the radial moments of the source 
function S. In spherical symmetry and with the line-element 
(14), the first two moment equations can be written (Thornc 
1981, equation [5.10c]) 



W2 — acso 



3a 3 & v 

+ — j(w2ar 3 )u — acsi = . 
br° 



(22a) 



(226) 



In equations (22) wo and wi have the dimensions of energy 
density and cso and csi are in units of erg cm~ 3 s _1 . It is 
well-known that the moment equations form a recursive sys- 
tem of differential equations that is not closed. At any given 
order k max it contains moments up to order k max +i in the 
frequency-integrated case (and k ma x+2 in the frequency- 
dependent case). This means that, in order to use these 
equations for calculations, it is necessary to make some "ad 
hoc" assumption to close the system (see e.g. Fu 1987, Cer- 
nohorsky & Bludman 1994 and references therein) and this 
is usually done on physical grounds by introducing suitable 
closure functions which relate Wfe mox+1 (and Wk max+2 where 
necessary) to moments of lower order. Since the behaviour of 
all of the moments is known in the asymptotic limits (when 
the interaction between matter and radiation is either very 
strong or completely absent), it is sufficient to prescribe a 
reasonable smooth function that connects these two limits 
(see e.g. NTZ). Clearly uncertainties in this will introduce 
some error into the calculation of the lower order moments, 
whose magnitude will be dependent on the closure relation 
but turns out to be no larger than ~ 15 % for the range of 
parameter values typical for real astrophysical flows (Turolla 
& Nobili 1988). 

The radiation hydrodynamic equations (15)-(21), to- 
gether with the first two moment equations (22) need to be 
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supplemented with the constitutive equations for the gas, 
the expressions for the source moments, a prescription for 
the closure function and suitable boundary conditions. 



A : 



(1.42 x l(T 27 T 1/2 /3 re! + 6.0 x i()- 22 T~ 1/2 ) 1 



h io 25 ( — "L — ^ 



3 -1 

erg cm s 



3 THE MODEL 

In the following, we will consider spherical accretion of a 
self-gravitating hydrogen gas in the gravitational field of a 
non-magnetized, non-rotating black hole. The basic equa- 
tions have been presented in the previous section; here we 
will specify all of the input physics (expressions for the 
source moments, equations of state, etc.) needed for solv- 
ing the problem. Stationary, spherical accretion onto black 
holes has recently been investigated in detail by NTZ. The 
main goal of the present paper is to ascertain the stability 
properties of the solutions found by NTZ; in particular, wc 
want to study the behaviour of the models in a certain range 
of accretion rates, for which both low and high luminosity 
solutions exist. To allow a direct comparison of our results 
with those of NTZ, we will adopt the same input physics 
which they considered, including the simplifying assump- 
tions that Compton scattering can always be treated in the 
Kompaneets limit and that pair processes can be neglected. 
This turns out to be an excellent approximation for LL mod- 
els while it becomes questionable in HL solutions where the 
gas temperature reaches ~ 10 10 K in the inner region. On 
the other hand, the work by Park (1990b) has shown that 
the inclusion of pair production and annihilation does not 
produce any dramatic changes. 

If the dominant radiative processes are bound-free, 
free-free and isotropic scattering, the radial source moments 
so and si can be written (Thorne 1981; NTZ) 



0r 



(1+4.4 x 10~ 10 T) 



so = p (e — fco»o) + kesPWo 



4k B T 
rn e c 2 



Si 



(23) 



(24) 



where ce is the emissivity per unit mass per unit time, ko, ki 
and k ea are the absorption, flux-mean and scattering opac- 
ities, respectively, and T 7 is the radiation temperature, de- 
fined by 



1 J °° hvwo{v)dv 
4k B f °° w (v)dv 



(25) 



In equation (23) the second term on the right hand side 
accounts for the energy exchange between matter and radi- 
ation due to non-conservative scatterings and is obtained by 
integrating the Kompaneets equation over frequency and ne- 
glecting the non-linear term which describes induced emis- 
sion. Since it is derived in the Fokker-Planck approximation, 
this term is certainly not adequate for describing the interac- 
tion between photons and electrons when the latter become 
relativistic, as happens in some of the solutions which we 
have computed. However, even in this case, the model can 
give a good qualitative indication of the correct results. We 
have expressed the emissivity using the interpolation of the 
cooling function A given by Stellingwerf & Buff (1982) 



P A 



(26) 



which includes bound-bound, free-bound, e—p and e-e 
brcmsstrahlung for a pure hydrogen gas; the factor /3 re ; is a 
relativistic correction. Assuming LTE between emitters and 
absorbers, we can use the Kirchhoff law to obtain the Planck 
mean opacity, kp = e/apT 4 ', where ap is the blackbody ra- 
diation constant. Since the actual spectral distribution of 
wo cannot be computed here, we use kp in place of the ab- 
sorption opacity fco- The flux-mean opacity k\ can be split 
into two terms: the first is the scattering opacity k es and the 
second is the sum of the contributions from all of the other 
radiative processes; however, since k es is always dominant 
for the range of densities and temperatures encountered in 
the present problem, we have approximated the additional 
term using the Rosseland mean kp calculated taking into 
account only free-free processes 



ki ~ k es + k R , 

kp = 6.4 x 10 22 pT~ 7/2 cm 2 g" 1 . 



(27) 



In the frequency-integrated transfer problem, the radiation 
temperature T 7 cannot be directly computed from its defini- 
tion (equation [25]). However, since T 7 appears only in the 
term in so which accounts for comptonization, it only be- 
comes important when the energy exchange between matter 
and radiation due to non-conservative scatterings starts to 
be effective. In a medium at rest, the fractional change of 
the mean photon energy (E = 4A;bT 7 ) because of scatterings 
with a thermal, non-relativistic distribution of electrons, fol- 
lows the relation 
dE 4k B T ( E 



Uk B T J 



1 adr . 



m e c-° \4k B T 

(Rybicki & Lightman 1979, Wandel, Yahil & Milgrom 1984) 
where r is the scattering depth, 4kpT(E /4kpT — l)/m e c 2 is 
the mean energy change per scattering and adr is the mean 
number of scatterings which a photon undergoes between r 
and t + dr, with 



a = 1 



a = 2r 



for r < 1 . 



for r > 1 . 



From the computational point of view, it is convenient to 
write an equation for T 7 which is valid over all of the integra- 
tion domain with continuous coefficients; for this reason the 
previous equation is usually written in the following form, 
which is an approximation near to r ~ 1 (Park & Ostriker 
1989; Park 1990a): 



dlnT 7 
din r 



2Y C 



(28) 



where Y c = 4feTmax(T,T 2 )/m e c 2 is the Compton param- 
eter. In stationary calculations, this equation has been used 
directly to give the variation of T 7 with r, at constant Eu- 
lerian time t, but for non-stationary flows it is not satis- 
factory to integrate it along the time-slice (i.e. at constant 
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Lagrangian time t) as this would imply an infinite speed of 
propagation of information. Instead, we apply it along the 
outward-pointing characteristics of the radiation field, fj, c (t) 
(defined from the moment equations [22]), and calculate the 
optical depth r along the same lines using 



f 



kespbFdfi . 



(29) 



This seems a reasonable choice because the radiation tem- 
perature is strictly related, by definition, to the radiation 
energy density and we expect that information will propa- 
gate along the characteristic lines of the radiation field. In 
this case it is not difficult to show (see the Appendix, equa- 
tion [A5]) that 



dt 



-acTvc 



(t- 1 ) 



i dr 7 



(30) 



where v c = bfi c /ac is the characteristic velocity for the ra- 
diation field and fi c — dfi c /dt. This is the actual form of 
the equation for T 7 used in our calculations. Equation (30) 
applies when comptonization is the dominant radiative pro- 
cess and also gives the correct behaviour (an outgoing radi- 
ation temperature wave) when non-conservative scattering 
becomes inefficient as long as true emission and absorption 
can be neglected. For the HL models, true absorption is 
never dominant and true emission is never likely to signif- 
icantly affect T 7 and so the use of equation (30) is always 
satisfactory. For the LL models, true emission and absorp- 
tion are dominant but in this case the second term on the 
right hand side of equation (23) is small compared with the 
other terms and an accurate evaluation of T 7 is no longer 
needed. The stationary limit of equation (30) (see again the 
Appendix, equation [A8]) is 



V VvJ dr 



(31) 



where the partial derivative is taken at constant Eulerian 
time t. In this form, the presence of a critical point where 
u = — Fv c is made apparent. We note that this result is a 
consequence of the finite velocity of propagation in equa- 
tion (30); in fact, as is well-known, the presence of critical 
points in the hydrodynamic equations for stationary flow 
is a relic of the characteristic velocity of the corresponding 
time-dependent equations. This result represents the main 
difference between the form of the T 7 equation used here and 
the one considered in all previous studies of this problem in 
the framework of black hole accretion (Park & Ostriker 1989; 
Park 1990a; NTZ). We stress that in all the discussion lead- 
ing to the equation for T 7 , we assumed a non-relativistic 
distribution for electrons. This has been done to ensure con- 
sistency with the Compton source term which has the simple 
expression given in equation (23) only in the non-relativistic 
limit. 

Finally, we need to specify the constitutive equations 
for a pure hydrogen gas 



P=[l+x(T)} 



pk B T 



(32) 



, = pc >\l + 3 -[x(T)+x*(T)]^ 



-[l-x(T)] 



E H 



(33) 



where T is the gas temperature, Eh = 13.6 eV is the hy- 
drogen ionization potential and x(T) is the degree of ioniza- 
tion, computed by equating the collisional-ionization and 
radiative-recombination rates (Buff & McCray 1974) and 
expressed using the interpolation formula of Stellingwerf and 
Buff (1982) 



x(T) 
with 



l + F 



(34a) 



exp 



In equation (33) 
x'(T) = l[0- 1 {v 



-1.58 x 10 5 K 



1) 



(346) 



where = k B T/rn e c 2 and i) = K 3 (0~ 1 )/K 2 {e~ 1 ) (K n is the 
modified Bessel function of order n). A polynomial fit by 
Service (1986) was used to calculate ij, giving an accuracy of 
a few parts in 10 5 . The third term inside the curly brackets in 
equation (33) accounts for the electrostatic potential energy 
of bound electrons in the neutral hydrogen atoms. 

The constitutive equations (32) and (33) can be used to 
express two fluid variables in terms of the other ones. Since 
the values of the temperature T and the density p are needed 
for evaluating the source moments so and si , it is more con- 
venient to use them in the hydrodynamic equations and to 
calculate e and P from equations (32) and (33). Substituting 
equation (33) into equation (15), the Energy equation can 
be written in the form 



+ p 



dx 
dlnT 

aso 



+ 



3 dx* 



2 dlnT 



k B T T t 
mpC 2 T 



(35) 



= 0. 



As far as the closure function is concerned, in the 
present calculation we chose to relate W2 to wo using the 
following expression 



/ = ^H[l+(-)T\ (36) 
w 3 L \t J J 

where to and n are free parameters. We made several trials 
with different expressions for / and found that the fractional 
difference between solutions obtained with different reason- 
able closures turns out to be no larger than ~ 20 %, which 
is acceptable here. In fact, a change in the closure param- 
eters was used to perturb the initial stationary solution, as 
we will discuss later. 



4 BOUNDARY CONDITIONS 

From the mathematical point of view, the equations of radi- 
ation hydrodynamics (16)-(21) and (35), the first two mo- 
ment equations (22) and the radiation temperature equation 
(30) form a hyperbolic system. In order for the problem to 
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Table 1. Parameters for the Stationary 
Models. 





m 


I 




r] = l/m 


1L 


28.5 


4.8 x 10- 


-6 


1.7 x 10- 


2L 


4.27 


9.2 x 10~ 


-7 


2.2 x 10" 


3L 


0.071 


3.5 x 10- 


-8 


4.9 x 10~ 


OH 


70.8 


7.3 x 10- 


-3 


1.0 x 10- 


1H 


28.3 


1.4 x 10- 


-3 


5.0 x 10- 


2H 


3.45 


2.2 x 10- 


4 


6.5 x 10- 



be well-posed, we need to specify values for all of the vari- 
ables at some initial time t — to over all of the integration 
domain pi n < p, < p ou t and also to assign suitable boundary 
conditions at the spatial boundaries pi n , p ut- The number 
of boundary conditions needed depends on the sign of the 
eigenvalues of the characteristic equations at each boundary. 
At the outer boundary, negative eigenvalues signify that in- 
formation is propagating into the integration domain from 
outside and a corresponding number of conditions must be 
assigned; the same is true for positive eigenvalues at the 
inner boundary. In the present case it can be shown that 
we need to prescribe 7 boundary conditions (4 at the inner 
boundary and 3 at the outer boundary) as follows: 2 condi- 
tions related to the fluid equations, 2 related to the moment 
equations, 1 each for the equations for T 7 , a and M. As far 
as the fluid boundary conditions are concerned: at /i — fiout 
we set a floating boundary (extrapolation in r) for u and, 
at p — fi in , we dropped the pressure gradient term from the 
Euler equation, making it advective in form and assuming 
strict free-fall very near to the black hole horizon. The in- 
ner conditions on T 7 and ti)i and the outer condition on wo 
were all taken as floating boundaries. The choice of a float- 
ing boundary is suitable when one does not want to put any 
constraint on a variable, leaving it free to adjust itself or to 
oscillate if there are waves propagating out of the integration 
domain (as for a vibrating string with free endpoints). 

As far as a and M are concerned, the time-slice at con- 
stant t is a characteristic direction for equations (19) and 
(20) and we put 



M = M 



at p = ii ou t , 



at fx — fiv 



(37) 
(38) 



The condition on a, equation (37), corresponds to synchro- 
nizing the coordinate time with the proper time of a comov- 
ing observer at the outer edge of the grid. This is also equal 
to the Eulerian time there, if the outer edge of the grid is 
placed sufficiently far away from the black hole. 

Finally, we note that, if the system tends to a stationary 
limit, the time-dependent equations reduce to their station- 
ary form and the solution which crosses any critical points 
in a regular way is automatically selected. 



5 NUMERICAL METHOD 

The equations of radiation hydrodynamics, presented in sec- 
tions 2 and 3, have been solved numerically for matter being 
accreted spherically onto a Schwarzschild black hole, using 
a Lagrangian finite difference scheme with a standard La- 
grangian organization of the grid. The code was adapted 
from one developed by Miller & Rezzolla (1995) for solution 
of the equations of radiation hydrodynamics in the context 



of the cosmological quark-hadron transition. We divided the 
integration domain (from the black hole horizon at r — r g 
out to 10 9 r g ) into a succession of comoving zones with each 
one having width Ap 21% larger than the one interior to 
it. Whenever the inner edge of the innermost zone crosses 
the horizon (which happens every 4-5 cycles with our time- 
step constraints) , we remove it from the calculation and per- 
form a regridding of all the variables. We followed a regrid- 
ding procedure previously adopted by Szuszkiewicz & Miller 
(1995) in connection with the study of disc accretion onto 
black holes. Originally a cubic spline interpolation was used, 
but this turned out to produce a numerically unstable evolu- 
tion in our case. A local cubic interpolation was eventually 
used instead, which was found to be satisfactory and effi- 
cient. At the initial time the effective mass contained within 
the inner boundary, Mo, is equal to the initial black hole 
mass Mbh, which we take to be 3M Q as in the stationary 
calculations. As time elapses, Mo increases as zones pass 
through the horizon and are removed from the calculation. 
However, during a characteristic evolutionary time interval, 
the mass of the material accreted is small compared with 
M bh . 

To have second-order accuracy in time, u and w\ are 
both evaluated at an intermediate time level. They are ad- 
vanced to the new time level at the end of each cycle after 
all of the other variables have been calculated. The time- 
step is adjusted in accordance with the relativistic Courant 
condition and two additional constraints on the fractional 
variations of p and T in each time-step, which are required 
to be smaller than 5%. In practice, we found that the time- 
step is usually limited by the last two conditions due to the 
fact that the variation of density and temperature, as seen 
by comoving observers, becomes very rapid near to the hori- 
zon where the flow velocity approaches the speed of light. 
As far as the spatial centering is concerned, p, T, wo, T 7 and 
a are treated as mid-zone quantities, while r, M, u and wi 
are treated as zone boundary quantities. 

Once the finite difference representation has been intro- 
duced, equations (16), (17), (19), (20) and (21) can be solved 
explicitly for u, p, a, M and r, respectively. Where necessary, 
linear interpolation and extrapolation in time were used to 
obtain the values of quantities at suitable time levels. The 
semi-logarithmic derivatives present in the Continuity equa- 
tion (17) and the equation for a (19) were solved using the 
Crank-Nicholson operator for equation (17) and the Leith- 
Hardy operator for equation (19) (see e.g. May & White 
1967). For the moment equations (22), we adopted a mixed 
representation: after dividing the first equation by wo and 
the second by mi we grouped together the terms in the fol- 
lowing way: 



{wp) t 
w 



( w i)t 



T-7» 



ac 

Wo 



■so 



/ 2 2x 



ac 

Wl 



Si 



:(fw ar 



(39a) 



(396) 



3ak {m>a4) »- abr^ 

where 102 has been expressed using the closure relation 
W2 = fwo with / being defined as in equation (36). The 
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Table 2. Time-scales. 

tj = — dynamical time 

uc - J 

t sg = ^-j sound crossing time for the gas v s = i ( ^ ) 

*sr = sound crossing time for the radiation v c = (/ + 3) ^ 

t c = e ~ P g +P cooling time C = pc + fceaico ^J ) 

*h = e ~ p ?, +P heating time _ff = pcu>o ( &o + fcea 4fcB "^">' ) 

*t/t = e ~\H C -cf thermal time 



terms on the left hand side of equations (39) were treated 
using the Crank-Nicholson operator, while the quantities 
appearing on the right hand side were calculated at the cor- 
rect time level by means of interpolation or extrapolation 
where necessary. Because the dependence on temperature in 
the Energy equation is rather sensitive, we adopted a semi- 
implicit scheme for equation (35) using a secant iteration 
method. The temperature at the new time level is calcu- 
lated iteratively starting from two initial estimates, based on 
the value at the previous time-step. Convergence is rapidly 
achieved in 4-5 iterations. Since so is in general very sensi- 
tive to the value of the temperature, the iteration was ex- 
tended to include also the zero-th moment equation (39a). 
The equation for T 7 (equation [30]) presents a particular nu- 
merical problem because of the delicate balance between T 
and T 7 which, following Bowers & Wilson (1991), we treated 
using a fully implicit differencing for (T — T 7 ) to achieve nu- 
merical stability. 

6 NUMERICAL RESULTS AND DISCUSSION 

We have calculated the time evolution of 6 models, start- 
ing from the stationary solutions listed in table 1 whose 
position in the (m, /) plane is shown in Fig. 1. Because of 
the different form of the T 7 equation used here, the present 
stationary solutions differ slightly from those of NTZ (see 
equation [31] and the following discussion). The definitions 
of relevant time-scales for the present discussion are listed 
in table 2. Along with the dynamical time-scale, ta (which is 
the characteristic time for an element of fluid with velocity 
u to travel a distance r), we have listed also the sound cross- 
ing times for the gas, t sg , and for the radiation, t sr (these 
are defined as the ratio of the radial length scale to the rel- 
evant characteristic velocity and represent the time needed 
for a "sound" wave in the gas or radiation fluid to travel 
a distance r). For determining whether stationary solutions 
are stable or not, we should in principle evolve each model 
for a time comparable with the relevant t ag or t 3r in order 
to allow information to travel from the inner regions to the 
outer ones. At 10 5 r 9 , t sg ~ 10 4 s and t sr ~ 1 s. To reach 
an evolutionary time t — 10 4 s would require prohibitively 
high computational times but in fact, as we shall see, all of 
the models evolve on a much shorter time-scale (typically of 
a few seconds) which is mainly determined by the thermal 
and radiative processes. The thermal balance is regulated by 
the cooling and heating time-scales, t c and th, which are the 
ratios of the internal energy of the gas to the cooling rate 
(C) and heating rate (H), respectively, and are both defined 
in table 2. We have also introduced the thermal time-scale, 
t t h, defined as the ratio of the internal energy of the gas to 



the net rate of energy input or output for the gas, \H — C\. 

For each model we started the numerical calculation 
from an initial perturbed solution which was calculated by 
changing the closure function. By varying to and n in equa- 
tion (36), we obtained a perturbation of the order of 10-20 % 
in the gas temperature and in the radiation moments. This 
way of setting the perturbation was not effective for model 
3L, which is optically thin everywhere, and so, only in this 
case, we decided to evolve the solution without any initial 
perturbation just to have an indication of the intrinsic sta- 
bility of these models. The solution after 14 seconds remains 
exactly in the initial stationary state. This is not surprising 
since, in this case, we know that cooling is very inefficient 
and so the result obtained by Moncrief (1980), who found 
that adiabatic flows are stable, was likely to apply here. In 
fact, optically thin solutions are not of great interest and we 
did not spend further time in the numerical analysis of the 
stability of these models. 

In Fig. 2 the results from the numerical calculation for 
model 1L are shown. This solution is representative of the 
behaviour of all optically thick LL models. As can be seen 
from the figure, the solution relaxes toward the stationary 
state (shown with a continuous line) on a time-scale of the 
order of t t h which, for this solution, is much shorter than 1 
second within 10 3 r g . This shows that these solutions are sta- 
ble, in agreement with the result obtained by Vitello (1984). 
The perturbation does not directly involve velocity and den- 
sity, which remain essentially equal to their initial values and 
the accretion rate also remains extremely constant. Radia- 
tion and gas pressure have negligible effect in these cases 
and matter is essentially in free-fall from the sonic radius 
(located at r ~ 10 9 r 9 ) down to the black hole horizon. Tem- 
perature and luminosity relax very quickly to their station- 
ary values. In the optically thick inner core, matter and ra- 
diation are in LTE and the luminosity is proportional to 
the local value of the temperature gradient; in the outer re- 
gion compressional heating balances free-bound cooling and 
the gas is essentially in radiative energy equilibrium at the 
hydrogen recombination temperature, T ~ IQ^K. 

For models on the high luminosity branch, the be- 
haviour of the time-dependent solutions is completely dif- 
ferent. In Figs. 3 and 4, we show the results of numerical 
calculations for models 2H and 1H, respectively. The mat- 
ter fluid is essentially in free-fall and the accretion rate is 
roughly constant. The interpretation of the temperature and 
luminosity profiles is less straightforward. They do not seem 
to converge at all toward the stationary solutions found by 
NTZ, but, after 70 seconds, they are still in essentially the 
initial perturbed state. We postpone discussion of these so- 
lutions to the end of this section. 
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Figure 2 




Figure 2. The Mach number u/Tv s , the accretion rate in Eddington units rh, the gas temperature T and the radiative luminosity 
I = Anr 2 cwi / 'L e are plotted versus log(r/r 9 ) for model 1L at different times. 



For larger m (model 1H), a thermal instability appears 
in the inner part of the flow after about 2-3 seconds, as can 
be seen from Fig. 4, and the temperature increases by al- 
most an order of magnitude. The cause of the onset of this 
instability will be discussed later. A few seconds after this, 
the velocity profile starts to deviate significantly from free- 
fall owing to the large drag exerted by the internal pressure 
gradients. A compression wave develops, whose front be- 
comes progressively steeper as it propagates outward and, 
after 8-10 seconds, a hydrodynamic shock forms at around 
10 3 -10 4 r 9 , where the Mach number passes through unity 



(small-dashed line in Fig. 4). Across the shock, the kinetic 
energy of the gas is dissipated into thermal energy and the 
density increases; matter starts to accumulate at the shock 
front and a corresponding decrease in p is seen in the inner 
region. Immediately behind the shock front, the gas accel- 
erates and free-fall is rapidly restored. We note that since 
the shock is very far away from the black hole horizon, the 
kinetic energy dissipated there is relatively small and the ra- 
diative luminosity does not increase significantly through it. 
As far as we know, this is the first time that a shock has been 
found in self-consistent solutions of black hole accretion (see 
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Figure 3 




Figure 3. The same as Fig. 2 for model 2H. 



Chang & Ostriker 1985 for a discussion of shock formation 
in stationary models). The large increase of I in the first 8- 
10 seconds (by more than one order of magnitude) is due 
to the enhancement in efficiency of free-free and Compton 
cooling caused by the increase in T. At later times, I starts 
to decrease because of the fall in density which leads to a 
decrease in efficiency of the cooling processes interior to the 
shock radius. The luminosity profile then has the typical be- 
haviour shown in Fig. 4 with the long-dashed line. This time 
variation is seen by a distant observer as a significant initial 
transient increase lasting ~ 8 s followed by a relatively rapid 
decrease at later times. Looking carefully at the Mach num- 



ber profile, it can be seen that the shock front (where the 
Mach number falls below one) is moving outward, at an ap- 
proximate speed of 10 8 cm/s ~ lCT 2 c. Hence this solution is 
definitely non-stationary as confirmed by the fact that the 
accretion rate is not constant and matter keeps accumulat- 
ing at the shock front. The evolution of model 1H (Fig. 4) 
was followed inserting a source of artificial viscosity, with the 
aim of getting a better treatment of the shock region. Follow- 
ing the standard prescription by von Neumann and Richt- 
myer (1950), a dissipative term Q oc p;_i/2 (w; — Ui-i) 2 , was 
inserted into the equations (here the subscripts indicate the 
locations on the finite difference grid at which each variable 
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Figure 4 




Figure 4. The same as Fig. 2 for model 1H. 



is evaluated; i represents a zone boundary and i — 1/2 rep- 
resents a mid-zone). However, since the flow is being com- 
pressed continuously from the sonic radius down to the black 
hole horizon, the amount of dissipation would be excessive, 
especially in the vicinity of r g , unless some modification is 
made to the standard procedure. In view of this, we decided 
to switch on the artificial viscosity only when the fractional 
variation of u across a grid zone a = 2(ui — Ui-i)/(ui+Ui-i) 
becomes larger than 30%. As the shock forms, a increases 
above 0.3 and the viscous term 

Qi-1/2 = k (l - Pi-1/2 («i - Ui-if <? (40) 



then starts to be progressively more effective. Here k is an 
adjustable constant (k — 2 in the actual calculation). About 
30 seconds after the beginning of the evolution and approx- 
imately 20 seconds after the shock formation, we were nev- 
ertheless forced to stop the calculation because of the for- 
mation of large numerical oscillations at the shock front and 
some more sophisticated treatment would clearly be desir- 
able. However, for the present purposes, we were content 
simply to demonstrate the existence of the shocked solu- 
tions. 

The evolution described for model 1H, with the appear- 
ance of a thermal instability and the formation of a hydro- 
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dynamic moving shock, is a common feature of all of the 
high luminosity models along the high rh part of the HL 
branch. We have made a systematic search for the point on 
the (rh, I) plane which marks the onset of the instability 
and the result is shown in Fig. 1, where all of the HL un- 
stable models are plotted as open circles. According to the 
analysis of Field (1965) and Stellingwerf (1982), the form of 
the free-free cooling function implies that the gas should be 
thermally unstable to isobaric short-wavelength perturba- 
tions so that, if the Compton heating rate exceeds the cool- 
ing rate at some radius, it will continue to do so until matter 
there has been heated to a temperature which is essentially 
equal to that of the radiation. In the present case, owing 
to the large value of the Compton parameter Y c , Compton 
cooling is equally as efficient as Compton heating and the 
analysis by Stellingwerf (1982) does not strictly apply. How- 
ever, as discussed by Cowie, Ostriker & Stark (1978), the 
instability is clearly due to the fact that the heating rate 
is greater than the cooling rate and, at the same time, the 
heating time is shorter than the dynamical time. In Figs. 5 
and 6, we have plotted the ratios of the heating time (t h ) 
to the cooling time (t c ) and of the thermal time (tth) to the 
dynamical time (t d ). These quantities are plotted against 
r/r g at different times for models 2H and 1H, respectively. 

As can be seen from Fig. 6, at the beginning th is slightly 
smaller than t c in the region around 10 3 -10 4 r s . There, heat- 
ing exceeds cooling and, since the flow cannot advect the 
excess energy efficiently (td — tth), a small perturbation is 
sufficient to make heating more effective and the onset of 
a thermal instability is unavoidable. Also the value of the 
thermal time-scale, Uh(W 3 r g ) ~ td(10 3 r 9 ) ~ 1 s, is consis- 
tent with the time-scale for the onset of the instability found 
numerically. The region of instability then moves to larger 
radii, as can be seen from Fig. 6, and so it is not surprising 
that the shock keeps moving outward. On the other hand, 
Fig. 5 shows that the ratio tth/td is significantly larger than 
unity for model 2H in the region where heating is more effec- 
tive than cooling and so the thermal instability is advected 
into the hole on a dynamical time-scale. In other words, 
since the models along the low rh part of the HL branch have 
smaller gas densities, the radiative heating and cooling are 
comparatively less efficient than compressional heating and 
the gas is essentially adiabatic. Only very far away, around 
3 x 10 5 r 9 , do the conditions seem to be favourable for the 
onset of the thermal instability; however, the thermal time- 
scale in that region is ~ 10 3 s and the evolutionary time 
becomes very large. This means that these solutions might 
also be unstable but on much longer time-scales. 

Finally, we note that, during the evolutionary times con- 
sidered, we did not find any evidence for possible transitions 
between the HL and the LL branch. 

7 CONCLUSIONS 

We have presented a systematic analysis of the stability and 
time-dependent behaviour of spherical accretion onto black 
holes in the framework of general-relativistic radiation hy- 
drodynamics. We have computed the evolution of a number 
of models from an intial perturbed state and confirmed that 
all of the low luminosity, NTZ stationary solutions, which 
are characterized by negligible comptonization, are indeed 
stable to thermal and radiative perturbations in agreement 



with previous investigations (e.g. Gilden & Wheeler 1980; 
Vitello 1984). On the other hand, the time evolution of 
high luminosity solutions, for which self-comptonization of 
bremsstrahlung photons is the main radiative and thermal 
process, exhibits a much richer phenomenology. We find that 
the upper part of the HL branch (for rh <; 10) enters the re- 
gion of the (rh, I) plane where preheating effects start to 
be important and this leads to the onset of a strong ther- 
mal instability giving rise to the formation of an outward- 
propagating hydrodynamic shock. These shocked solutions 
show significant transient increases in the total luminosity. 
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Figure 5 




Figure 5. Ratios of the heating time (t/j) to the cooling time (i c ) and of the thermal time (tth) to the dynamical time (t^) plotted 
against log(r/r 9 ) for model 2H at different times. 
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Figure 6. The same as Fig. 5 for model 1H. Note that the curve for t = in the left panel is difficult to see because it is almost 
completely hidden by the continuous line. 



Rezzolla, L., & Miller, J.C. 1994, Class. Quantum Grav., 11, 1815 
Rybicki, G.B., &: Lightman, A. P. 1979, Radiative Processes in 

Astrophysics (New York: Wiley) 
Schmid-Burgk, J. 1978, Ap&SS, 56, 191 
Service, A.T. 1986, ApJ, 307, 60 
Shapiro, S.L. 1973a, ApJ, 180, 531 
Shapiro, S.L. 1973b, ApJ, 185, 69 
Shull, J.M. 1979, ApJ, 229, 1092 
Shvartsman, V.F. 1971, Soviet Astr.-AJ, 15, 377 
Soffcl, M.H. 1982, A&A, 116, 111 
Stellingwerf, R.F., & Buff, J. 1978, ApJ, 221, 661 
Stellingwerf, R.F., & Buff, J. 1982, ApJ, 260, 755 
Stellingwerf, R.F. 1982, ApJ, 260, 768 
Szuszkicwicz, E., & Miller J.C. 1995, in preparation 
Tamazawa, S., Toyama, K., Kancko, N., & 6no, Y. 1974, Ap&SS, 

32, 403 



Thorne, K.S. 1981, MNRAS, 194, 439 

Turolla, R., & Nobili, L. 1988, MNRAS, 235, 1273 

Vitello, P.A.J. 1978, ApJ, 225, 694 

Vitello, P.A.J. 1984, ApJ, 284, 394 

von Neumann, J., & Richtmyer, R.D. 1950, J. Appl. Phys., 21, 
232 

Wandel, A., Yahil, A., & Milgrom, M. 1984, ApJ, 282, 53 



APPENDIX A 

As discussed in the text, we take equation (28) to be valid 
along the outward-pointing characteristics of the radiation 
field, /J, c (t), with the optical depth being defined along the 
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same lines by 

k es pbYdp . 



-f 



(Al) 



(A2) 



Then, we have 

^- = -k ea pbF. 
dp 

(Throughout this Appendix the total derivatives are taken 
along the outward-pointing characteristics of the radiation 
field). Using equations (A2) and (28), we can write 

dT-, dT-,dr. , ,_, . dT-, 



dr 
1 



(A3) 



where p c = dp c /dt and v c = bp c /ac = (/ + 1/3) 1 / 2 . How- 
ever, also 



dT-, 
~df 



dT-, 

dt 



dT-, 
dp 



Pc 



and so we finally get 



dT-, 
dt 



-acTvc 



7 V t V + br dp 



(A4) 



.(A5) 



In practice, we obviously cannot calculate the integrated 
value of t directly from expression (Al) evaluated along the 
outward-pointing characteristics for the radiation since this 
would involve knowledge of information ahead of the current 
time reached in the calculation. Instead, we evaluated equa- 
tion (Al) along the time-slice and this should give reason- 
able values. While it is important to calculate the derivative 
(A2) along the correct directions in order to ensure a satis- 
factorily causal propagation of information, the calculation 
of the integral (Al) should not be so sensitive. 

To derive the stationary limit of equation ( A5) , we need 
first to write it in terms of the Eulerian time and radial 
coordinates (t,r). Using the chain rule for differentiation in 
equation (A5), we obtain 



tt- 



dT-, 
dt 



+ rf 



dT-, 
dr 



= — acTvc 



+ ±t^ 



1 dT-, 

+ &r r " dr 



(A6) 



The condition of stationarity is expressed with respect to 
the fixed (constant r) Eulerian observer by 



d_ 

dt 



0. 



(A7) 



Taking the stationary limit in equation (A6) and using r t = 
acu and r M = bT finally yields 



dT-, 
dr 



^kcspYc ^ / T 



(!->) 



(A8) 



which corresponds to equation (31). 
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